Sizing-up the WIMPs of Milky Way : Deriving the velocity distribution of Galactic Dark 

Matter particles from the rotation curve data 



Pijushpani Bhattacharjee 1 *, Soumini Chaudhury 1 ^, Susmita Kundu 1 "' 1 and Subhabrata Majumdar 2 ^ 
AstroParticle Physics & Cosmology Division and Centre for AstroParticle Physics, 
Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700064- India 
2 Department of Theoretical Physics, Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005. India 

The velocity distribution function (VDF) of the hypothetical Weakly Interacting Massive Particles 
(WIMPs), currently the most favored candidate for the Dark Matter (DM) in the Galaxy, is deter- 
mined directly from the rotation curve data of the Galaxy assuming isotropic VDF. This is done by 
"inverting" — using Eddington's method — the Navarro-Frenk- White universal density profile of 
the DM halo of the Galaxy, the parameters of which are determined, by using Markov Chain Monte 
Carlo (MCMC) technique, from a recently compiled set of observational data on the Galaxy's rota- 
tion curve extended to distances well beyond the visible edge of the disk of the Galaxy. The derived 
most-likely local isotropic VDF strongly differs from the Maxwellian form assumed in the "Stan- 
dard Halo Model" (SHM) customarily used in the analysis of the results of WIMP direct-detection 
experiments. A parametrized (non-Maxwellian) form of the derived most-likely local VDF is given. 
The astrophysical "g-factor" that determines the effect of the WIMP VDF on the expected event 
rate in a direct-detection experiment can be lower for the most-likely VDF than that for the closest 
Maxwellian VDF by as much two orders of magnitude at the lowest WIMP mass threshold of a 
typical experiment. 



Several experiments worldwide are currently trying to di- 
rectly detect the hypothetical Weakly Interacting Massive 
Particles (WIMPs), thought to constitute the Dark Mat- 
ter (DM) halo of our Galaxy, by looking for nuclear recoil 
events due to scattering of WIMPs off nuclei of suitably cho- 

| sen detector materials in low background underground fa- 
cilities. The rate of nuclear recoil events depends crucially 

[ on the local (i.e., solar neighbourhood) density and veloc- 
ity distribution of the WIMPs in the Galaxy [1], which are 
a priori unknown. Estimates based on a variety of obser- 
vational data typically yield values for the local density of 
DM, p DM , Q , in the range 0.2 - 0.4 GeVcm" 3 (5.27 x 10~ 3 - 
0.01 Mq pc~ 3 ) [2]. In contrast, not much knowledge directly 
based on observational data is available on the likely form 
of the velocity distribution function (VDF) of the WIMPs 
in the Galaxy. The standard practice is to use what is of- 
ten referred to as the "Standard Halo Model" (SHM), in 
which the DM halo of the Galaxy is described as a single- 
component isothermal sphere [3], for which the VDF is as- 
sumed to be isotropic and of Maxwell-Boltzmann (hereafter 
simply "Maxwellian") form, /(v) cx cxp(— |v| / vo 2 ), with a 
truncation at an assumed value of the local escape speed, 
and with vq — u c .©, the circular rotation velocity at the 
location of the Sun. Apart from several theoretical issues 
(see, e.g., [4]) concerning the self-consistency of the SHM as 
a model of a finite-size, finite-mass DM halo of the Galaxy, 
high resolution cosmological simulations of DM halos [5] give 
strong indications of significant departure of the VDF from 
the Maxwellian. On the other hand, these cosmological sim- 
ulations do not yet satisfactorily include the gravitational 
effects of the visible matter components of the real Galaxy, 
namely, the central bulge and the disk, which provide the 
dominant gravitational potential in the inner regions of the 
Galaxy including the solar neighborhood region. 

The VDF of the DM particles at any location in the 
Galaxy is sclf-consistently related to their spatial density as 
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well as to the total gravitational potential, $(x), at that 
location. For a spherical system of collisionless particles 
(WIMPs, for example) with isotropic VDF satisfying the col- 
lisionless Boltzmann equation, the Jeans theorem [3] ensures 
that the phase space distribution function (PSDF), J-"(x, v), 
depends on the phase space coordinates (x, v) only through 
the total energy (per unit mass), E = ^v 2 + $(r), where 
v = |v|, r = |x|. For such a system, given a isotropic spatial 
density distribution p(r) = j d 3 vJ-(E), one can get a unique 
T by the Eddington formula [3, 6] 
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where \£(r) = — <j>(r) + $(r = oo) is the relative potential and 
£ = —E+$(r = oo) = ^(r)— \v 2 is the relative energy, with 
T > for £ > 0, and T — for £ < 0. The latter condition 
implies that at any location r, the VDF, / r (v) = T j p{r) , 
has a natural truncation at a maximum value of v, namely, 
iwO) = y/2^{r). 

Thus, given a isotropic density profile of a set of colli- 
sionless particles, we can calculate the VDF, / r (v), using 
equation (1) provided the total gravitational potential $(r) 
in which the particles move is known. A direct observational 
probe of $(r) is provided by the rotation curve (RC) of the 
Galaxy, the circular velocity of a test particle as a function 
of the galactocentric distance. In this paper we reconstruct 
the total gravitational potential $(r) in the Galaxy directly 
from the Galactic RC data and then use equation (1) to ob- 
tain the VDF, /r(v), of the WIMPs at any location in the 
Galaxy [7]. 

We shall assume that the DM density profile to be used 
on the right hand side of equation (1) is of the universal 
NFW [9] form, which, when normalized to DM density at 
solar location, p DM ,o, can be written as 
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where Rq is the distance of Sun from the Galactic centre. 
The profile (2) has two free parameters, namely, the density 
Pdm.o and the scale radius r s . 



The total gravitational potential seen by the DM particle, 
<I>, is given by $ = $ DM + <J> VM , where $ DM is the DM poten- 
tial corresponding to the density distribution (2) and 4> VM 
is the total potential due to the visible matter (VM) com- 
ponent of the Galaxy. The latter can be effectively mod- 
eled [10] in terms of a spheroidal bulge superposed on an 
axisymmetric disk, with density distributions given, respec- 
tively, by Bulge : p b = p b0 (l + (r/ r b ) 2 ) 3/2 , where p b0 and 
r b are the central density and scale radius of the bulge, re- 
spectively, and Disk : p d (R, z) = £si e -(.R- R@)/ R* e -M/*d ; 
where R and z are the axisymmetric cylindrical coordinates 
with r = (R 2 + z 2 ) 1 / 2 , R d and z d are the scale length and 
scale height of the disk, respectively, and E© is its local sur- 
face density. The corresponding gravitational potentials for 
these density models, $buige and $disk, can be easily obtained 
by numerically solving the respective Poisson equations, giv- 
ing $ VM = $bulgc + *disk- 

The density models specified above have a total of seven 
free parameters, namely, r s , p OM ,@, p b «, r h , E Q , i? d , and 
z d . We determine the most-likely values and the 68% C.L. 
upper and lower ranges of these parameters by performing 
a Markov Chain Monte Carlo (MCMC) analysis (see, e.g., 
Refs. [11]) using the observed RC data of the Galaxy. For 
a given set of the Galactic model parameters, the circular 
rotation speed, v c (R), as a function of the Galactocentric 
distance R, is given by 



v 2 c (R) = R-^ [* DM 0R, z = 0) + $ VM (i?, z = 0) 



(3) 



For the observational data, we use a recently compiled set 
of rotation curve data [12] that extends to Galactocentric 
distances well beyond the visible edge of the Galaxy. This 
data set corresponds to a choice of the Local Standard of 
Rest (LSR) set to (R e ,v Ct& ) = (8.0kpc, 200kms~ 1 ) [13]. 
For the MCMC analysis, we use the x 2 -test statistic defined 

as X 2 = ( " C '^ S c„!r th J ' Where ^,obs and <crror a re, 

respectively, the observational value of the circular rotation 
speed and its error at the i-th value of the galactocentric dis- 
tance, and v\ th is the corresponding theoretically calculated 
rotation speed. For priors on the free parameters involved, 
we have taken the following ranges of the relevant parameters 
based on currently available observational knowledge : For 
the VM parameters, p b0 : [0.1 - 2] x 4.2 x 10 2 M© pc -3 [10]; 
r b : [0.01-0.2] x 0.103 kpc [10]; E© : [35 — 58] M© pc~ 2 [14]; 
R d : [1-7— 3.5] kpc [10, 15]. The parameter z d has been 
fixed at 340 pc [16] since the results are fairly insensitive to 
this parameter. For the DM parameters we took a wide 
enough prior range for r s : [0.1 — 100] kpc and /9 D m,© : 
[0.1 — 0.5] GeVcm~ 3 consistent with values recently quoted 
in literature [2] . 

The results of our MCMC analysis are summarized in Ta- 
ble I and Figure 1. Figure 2 shows the theoretically cal- 
culated rotation curve for the most-likely set of values of 
the Galactic model parameters obtained from the MCMC 
analysis and listed in Table I, and its comparison with the 
observed rotation curve data. In Table II, we display the 
values of some of the physical quantities of interest charac- 
terizing the Galaxy, derived from the Galactic parameters 
listed in Table I. The values in Table II are in reasonably 
good agreement with the values of these quantities quoted 
in recent literature [8, 12, 17]. The relatively large uncer- 
tainties in the values of some of the quantities that receive 
dominant contribution from the DM halo properties at large 
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0.092 


54.30 


3.14 


SD 
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TABLE I: The most-likely values of the Galactic model param- 
eters, as well as their 68% C.L. lower and upper ranges, means 
and standard deviations (SD), obtained from our MCMC analysis 
using the observed rotation curve data. 
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FIG. 1: The 2D posterior probability density function for Dark 
Matter parameters (r a — Pdm,q), marginalized over the visible 
matter parameters. 



Galactocentric distances are simply a reflection of the rela- 
tively large uncertainties of the rotation curve data at those 
distances. 

The Galactic model parameters determined above allow 
us to reconstruct the total gravitational potential $(x) at 
any location in the Galaxy. Because of the axisymmetric 
nature of the VM disk, this potential is non-spherical. To 
use equation (1), which is valid only for a spherical sym- 
metric situation, we use the spherical approximation [8, 18], 
$ VM (r) ~ G J r M VM (r')/r' 2 dr', where M V m is the total VM 
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FIG. 2: Rotation curve of the Galaxy with the most-likely 
values of the Galactic model parameters listed in Table I 
data with error bars are from Ref. [12]. 
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TABLE II: The most-likely values of various relevant physical pa- 
rameters of the Milky Way and their upper and lower ranges de- 
rived from the most-likely- and 68% C.L. upper and lower ranges 
of values of the Galactic model parameters listed in Table I. 



mass contained within r. 

The resulting normalized speed distribution, f r (v) = 
(inv 2 / p DM (r)) / r (v) (with J f r (v)dv = 1) evaluated at the 
location of the Sun, giving fo(v), is shown in Figure 3. For 
comparison, we also show in the same Figure the results from 
various large N-body simulations as well as the distribution 
corresponding to the closest Maxwcllian form, axwe ^ (v) oc 
v 2 exp (— v 2 / Vq) , truncated at Umax,© = SlGkms^ 1 (see 
Table II), with the free parameter vq determined to be 
206 kms- 1 . 

As evident from Figure 3, the speed distribution differs 
significantly from the Maxwellian form. We find that the 
following parametrized form, which goes over to the standard 
Maxwellian form in the limit of the parameter k — > 0, gives a 
good fit to our numerically obtained most-likcly local speed 
distribution shown in Figure 3: 



f Q (v) « (4™ 2 /Pdm, ) (£(/3) - £(/3 ma *)) , 
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where ^(x) 
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(1 + x) k e"* , p 
= 339kms _1 and k = 



-- v 2 /v 2 , /3 max = 
T.47. As a quan- 



u max,0 

titativc measure of the deviation of a model form of the lo- 
cal speed distribution, j™ 0010 ^ from the numerically obtained 
most-likely (ML) form, / ML , shown in Figure 3, the quan- 

tity X j = (W)Eti [/ ML (^) - / modcl (^)] 2 has a value 
of ~ 7.2 x 10~ 5 for the parametrized form (4) compared to 
a value ~ 1.7 x 10 -3 for the closest Maxwellian shown in 
Figure 3. Note also that our results differ significantly from 
those obtained from the N-body simulations. 

In Figure 4 we show the most-likely f r (v)'s at several dif- 
ferent values of the Galactocentric distance r. Notice how 
the peak of the distribution shifts towards smaller values 
of v and the width of the distribution shrinks, as we go to 
larger r, with the distribution eventually becoming a delta 
function at zero speed at asymptotically large distances, as 
expected. The non-Maxwellian nature of the distribution at 
all locations is also clearly seen, with the Maxwellian ap- 
proximation always overestimating the number of particles 
at both low as well as extreme high velocities. The inset in 
Figure 4 shows our results for the pseudo phase space den- 
sity, Q = p/(v 2 ) 3 / 2 , as a function of r, and its comparison 
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FIG. 3: Normalized local speed distribution, /©(w), correspond- 
ing to the most-likely set of values of the Galactic model pa- 
rameters given in Table I (black curve) and its uncertainty band 
(yellow) corresponding to the 68% C.L. upper and lower ranges 
of the Galactic model parameters. The magenta curve (almost 
overlapping with the black one) corresponds to the parametrized 
fit to the most-likely distribution given in equation (4). The 
closest Maxwellian, /g axwoll (?;) oc v 2 exp (-v 2 /vq), truncated at 
■Umax,© = 516 km s -1 (see Table II), with the free parameter vo 
determined to be 206 km s -1 , as well as results from some of the 
large N-body simulations [5], are also shown for comparison. 
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FIG. 4: Normalized speed distribution of the DM particles at 
various Galactocentric radii (solid curves), corresponding to the 
most-likely set of values of the Galactic model parameters given 
in Table I. The curves (dotted) for the corresponding closest 
Maxwellians are also shown for comparison. The inset shows the 
pseudo-phase space density of DM, Q = p/(v 2 ) 3 ^ 2 , as a function 
of r. 



with the power-law behavior predicted from simulation re- 
sults [19]. Note the agreement with the power-law behavior 
at large distances but strong deviation from it at smaller 
Galactocentric radii, which we attribute to the effect of the 
visible matter: For a given DM density profile, the additional 
gravitational potential provided by the VM supports higher 
velocity dispersion of the DM particles, making Q smaller 
than that for the DM-only case. 

We now discuss the implications of our results for the anal- 
ysis of direct detection experiments. The differential rate of 
nuclear recoil events per unit detector mass (typically mea- 
sured in counts/day/kg/keV), in which a WIMP (hereafter 
generically denoted by x with mass m x ) elastically scatters 



4 



off a target nucleus of mass mjy leaving the recoiling nucleus 
with a kinetic energy E R , can be written as [1] 



lowest WIMP mass that can 
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p x g(E R ,t), (5) 



where p x = /o DM ,o is the mass density of WIMPs in the 
solar neighborhood, a(q 2 ) is the momentum transfer de- 
pendent effective WIMP-nucleus elastic cross section, fi = 
m x mN I (m x +mN) is the reduced mass of the WIMP-nucleus 
system, and 



g(E R ,t) 



v E (t)) 0(u max - u min ) , 
(6) 

is the crucial "g-factor" that contains all information about 
the local VDF of the WIMPs [20]. In (6) the variable u 
(with u = |u|) represents the relative velocity of the WIMP 
with respect to the detector at rest on Earth, and vs(t) is the 
(time-dependent) velocity of the Earth relative to the Galac- 

1 /2 

tic rest frame. The quantity u m i n (E R ) = (rriNE R /2fj, 2 ) is 
the minimum WIMP speed required for giving a recoil en- 
ergy E R to the nucleus, and u max (i) is the (time-dependent) 
maximum WIMP speed [4] corresponding to the maximum 
speed f m ax (defined in the Galactic rest frame) for the VDF 
under consideration. 




FIG. 5: The ratio (solid curves), ( = gMh(E t h) / ffMaxwei^-Eth) , of 
the g-factor calculated with our most-likely (ML) form of /o(i>) 
shown in Figure 3 to that for the closest Maxwellian form also 
shown in Figure 3, as a function of the WIMP mass m x , for two 
different target nuclei, namely, Sodium and Xenon, both with 
E t h = 2keV. The shaded bands correspond to the uncertainty 
bands of /o(i>) shown in Figure 3. The calculations are for 2nd 
June, when the Earth's velocity in the Galactic rest frame is max- 



Note that the quantity g(E R ,t) takes its largest value at 
E R = E t h, the threshold energy for the experiment under 
consideration. To illustrate the effect of the non-Maxwellian 
nature of the VDF and its uncertainty, we define the quan- 
tity C = SMiX-EthV^MaxwciK-Eth), the ratio of the g-factor 
calculated with our most-likely (ML) form of /©(f) shown in 
Figure 3 to that for the closest Maxwellian form also shown 
in Figure 3, both evaluated at E R = E t h- A plot of C as a 
function of the WIMP mass m x , for two different target nu- 
clei, viz. Sodium and Xenon, in both case with E^ = 2 keV, 
is shown in Figure 5. 



experiment is given by 
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Figure 5, the effect of the departure from Maxwellian 
distribution is most significant at the lowest WIMP mass 
where the difference can be as much as two orders of 
magnitude. 



To summarize, a first attempt has been made to derive the 
velocity distribution (assumed isotropic) of the dark matter 
particles in the Galaxy directly using the rotation curve data. 
The distribution is found to be significantly non-Maxwellian 
in nature, the implication of which is a sizable deviation of 
the expected direct detection event rates from those calcu- 
lated with the usual Maxwellian form. 
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